Abundancia Relativa
Instalar y cargar paqueterías
if(!requireNamespace("tidyverse", quietly=TRUE))
install.packages("tidyverse")
library(tidyverse)
library(plotly)Importar los datos a analizar
De las tablas de abundancia de los distintos amplicones, generadas con el procesamiento de dada2, se importan los archivos generados.
data <- read.delim("~/Downloads/secuencias/Tablas/ASVs_total_silva.csv") #Adecuar la ruta del archivoObtener abundancia relativa de NA
Para tener una idea del porcentaje de lecturas que van a ser modificadas debido a la falta de categorización taxonómica, se obtiene la abundancia de NAs a lo largo de todas las muestras y ASVs.
# Sumar el conteo de NA y el total de lecturas para todas las muestras
total_NA <- data %>%
filter(is.na(Family)) %>%
select(starts_with("AC")) %>%
pivot_longer(everything(), names_to = "Sample", values_to = "Count") %>%
summarize(Total_NA = sum(Count))
# Sumar el total de lecturas en todas las muestras
total_reads <- data %>%
select(starts_with("AC")) %>%
pivot_longer(everything(), names_to = "Sample", values_to = "Count") %>%
summarize(Total_Reads = sum(Count))
# Calcular la abundancia relativa de NA en todo el conjunto de muestras
abundancia_relativa_NA <- total_NA$Total_NA / total_reads$Total_Reads * 100
# Mostrar el resultado
abundancia_relativa_NACalcular abundancia relativa
Familia
Para calcular la abundancia relativa, primero se deben reemplazar los valores nulos NAs por la categoría “No Clasificados”. Posteriormente se calcula la abundancia relativa por familia taxonómica.
data$Family <- ifelse(is.na(data$Family), "No Clasificados", as.character(data$Family))
family_abundance <- data %>%
select(starts_with("AC"), Family) %>%
pivot_longer(cols = -Family, names_to = "Sample", values_to = "Count") %>%
group_by(Sample, Family) %>%
summarize(Total = sum(Count)) %>%
mutate(Percentage = Total / sum(Total) * 100)A los ASVs con una abundancia menor al 1% se les asigna en la nueva categoría de Otros.
family_abundance <- family_abundance %>%
mutate(Family = ifelse(Percentage < 1, "Otros", Family))Posteriormente se procede a realizar de nuevo el cálculo de abundancia relativa.
family_abundance <- family_abundance %>%
group_by(Sample, Family) %>%
summarize(Total = sum(Total)) %>%
mutate(Percentage = Total / sum(Total) * 100)Género
Para los taxones de Género y Especie se sigue un paso adicional. Se reemplazan los valores nulos NA con el nombre del taxón Familia usando el formato “u_Nombre”.
##Se usa de nuevo la tabla original sin la categoría de Otros para no perder esas asignaciones taxonómicas
data <- read.delim("~/Downloads/secuencias/Tablas/ASVs_total_silva.csv")
## Reemplazar los valores NA en la columna "Genus" usando el taxon anterior.
data$Genus <- ifelse(is.na(data$Genus), paste0("u_", data$Family), data$Genus)
## Reemplazar "u_NA" y "u_Otros" con NA en la columna "Genus"
data$Genus <- ifelse(data$Genus == "u_NA" | data$Genus == "u_Otros", NA, data$Genus)
# Reemplazar los valores NA en la columna "Genus" con "No Clasificados"
data <- data %>%
mutate(Genus = if_else(is.na(Genus), "No Clasificados", Genus))genus_abundance <- data %>%
select(starts_with("AC"), Genus) %>%
pivot_longer(cols = -Genus, names_to = "Sample", values_to = "Count") %>%
group_by(Sample, Genus) %>%
summarize(Total = sum(Count)) %>%
mutate(Percentage = Total / sum(Total) * 100)genus_abundance <- genus_abundance %>%
mutate(Genus = ifelse(Percentage < 1, "Otros", Genus))genus_abundance <- genus_abundance %>%
group_by(Sample, Genus) %>%
summarize(Total = sum(Total)) %>%
mutate(Percentage = Total / sum(Total) * 100)Especies
Por último se puede realizar lo mismo con el taxón de especies para mayor especificidad.
##Se usa de nuevo la tabla original sin la categoría de Otros para no perder esas asignaciones taxonómicas
data <- read.delim("~/Downloads/secuencias/Tablas/ASVs_total_silva.csv")
## Reemplazar los valores NA en la columna "Species" usando el taxon de familia.
data$Species <- ifelse(is.na(data$Species), paste0("u_", data$Family), data$Species)
## Reemplazar "u_NA" y "u_Otros" con NA en la columna "Species"
data$Species <- ifelse(data$Species == "u_NA" | data$Species == "u_Otros", NA, data$Species)
# Reemplazar los valores NA en la columna "Species" con "No Clasificados"
data <- data %>%
mutate(Species = if_else(is.na(Species), "No Clasificados", Species))species_abundance <- data %>%
select(starts_with("AC"), Species) %>%
pivot_longer(cols = -Species, names_to = "Sample", values_to = "Count") %>%
group_by(Sample, Species) %>%
summarize(Total = sum(Count)) %>%
mutate(Percentage = Total / sum(Total) * 100)species_abundance <- species_abundance %>%
mutate(Species = ifelse(Percentage < 1, "Otros", Species))species_abundance <- species_abundance %>%
group_by(Sample, Species) %>%
summarize(Total = sum(Total)) %>%
mutate(Percentage = Total / sum(Total) * 100)Graficar abundancias
Para poder visualizar las abundancias relativas, se procede a graficar los diferentes taxones. En las gráficas generadas, Otros encapsula los ASVs con taxonomía no encontrada en la base de entrenamiento, y aquellos cuya abundancia relativa era menor al 1%.
Familia
family_abundance$Family <- factor(family_abundance$Family,
levels = c(setdiff(unique(family_abundance$Family), "Otros"), "Otros"))
family_plot<-ggplot(family_abundance, aes(x = Sample, y = Percentage, fill = Family)) +
geom_bar(stat = "identity", position = "stack") +
labs(title = "Abundancia Relativa por Familia", y = "Porcentaje", x = "Muestra") +
theme(axis.text.x = element_text(angle = 90, hjust = 1),
legend.key.size = unit(0.5, "cm"),
legend.text = element_text(size = 8),
plot.margin = margin(1, 1, 1, 1, "cm")
) + scale_color_brewer(palette="Set3") + guides(fill = guide_legend(ncol = 2))
ggplotly(family_plot)Género
genus_abundance$Genus <- factor(genus_abundance$Genus,
levels = c(setdiff(unique(genus_abundance$Genus), "Otros"), "Otros"))
genus_plot<-ggplot(genus_abundance, aes(x = Sample, y = Percentage, fill = Genus)) +
geom_bar(stat = "identity", position = "stack") +
labs(title = "Abundancia Relativa por Género", y = "Porcentaje", x = "Muestra") +
theme(axis.text.x = element_text(angle = 90, hjust = 1),
legend.key.size = unit(0.5, "cm"),
legend.text = element_text(size = 8),
plot.margin = margin(1, 1, 1, 1, "cm")
) +scale_color_brewer(palette="Set3") + guides(fill = guide_legend(ncol = 2))
ggplotly(genus_plot)Especies
species_abundance$Species <- factor(species_abundance$Species,
levels = c(setdiff(unique(species_abundance$Species), "Otros"), "Otros"))
species_plot<-ggplot(species_abundance, aes(x = Sample, y = Percentage, fill = Species)) +
geom_bar(stat = "identity", position = "stack") +
labs(title = "Abundancia Relativa por Especies", y = "Porcentaje", x = "Muestra") +
theme(axis.text.x = element_text(angle = 90, hjust = 1),
legend.key.size = unit(0.5, "cm"),
legend.text = element_text(size = 8),
plot.margin = margin(1, 1, 1, 1, "cm")
) + scale_color_brewer(palette="Set3") + guides(fill = guide_legend(ncol = 2))
ggplotly(species_plot)Guardar gráficas
Para guardar las gráficas en distintos formatos, se tiene que asignar una variable a la función ggplot. A continuación se muestra un ejemplo de como quedarían las líneas de código para realizar esto.
path<-"/Users/adriantorres/Downloads/secuencias"
ggsave(file.path(path, "family_abundance.svg"), plot = family_plot, dpi = 600)
ggsave(file.path(path,"family_abundance.jpg"), plot = family_plot, dpi = 600)
ggsave(file.path(path,"family_abundance.pdf"), plot = family_plot, dpi = 600)
#En este caso family_plot corresponde al nombre de la variable asignada a la función ggplot. Se realiza lo mismo con el taxón género y especies.
ggsave(file.path(path, "genus_abundance.svg"), plot = genus_plot, dpi = 600)
ggsave(file.path(path,"genus_abundance.jpg"), plot = genus_plot, dpi = 600)
ggsave(file.path(path,"genus_abundance.pdf"), plot = genus_plot, dpi = 600)ggsave(file.path(path, "species_abundance.svg"), plot=species_plot, dpi = 600)
ggsave(file.path(path,"species_abundance.jpg"), plot=species_plot, dpi = 600)
ggsave(file.path(path,"species_abundance.pdf"), plot=species_plot, dpi = 600)